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Abstract 



We present algorithms to evaluate two types of multiple sums, which appear in 
higher-order loop computations. We consider expansions of a generalized hyper- 



e-eometrir-tvne sums V T{ai-n+c 1 )T{a2-n+c 2 )---T{a M -n+c M ) m . . . n N 

geometric type sums, 2^n u -,n N r(b 1 - n +d 1 )r(b2-n+d2)---r(b M - n +d M ) x i X N 



Type I sum corresponds to the case where, in the limit e — > 0, the summand 
reduces to a rational function of n/s times x™ 1 • • -x^f; Cj,dj's can depend on an 
external integer index. Type II sum is a double sum (iV = 2), where Ci,di's are 
half- integers or integers as e — > and Xi = 1; we consider some specific cases 
where at most six T functions remain in the limit e — > 0. The algorithms enable 
evaluations of arbitrary expansion coefficients in e in terms of Z-sums and multiple 
polylogarithms (generalized multiple zeta values). We also present applications of 
these algorithms. In particular, Type I sums can be used to generate a new class 
of relations among generalized multiple zeta values. We provide a Mathematica 
package, in which these algorithms are implemented. 
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1 Introduction 



During the last decades, there have been remarkable developments in technologies for 
computing higher-order radiative corrections in quantum field theories. Applications 
of these computational technologies extend from phenomenology of particle physics to 
various fields of theoretical research; see, for instance, reviews [TJ [21 [3J, 111 El [7] . 

It has become a standard procedure to reduce numerous integrals, which appear in 
higher-order corrections, to a small set of integrals (master integrals) using identities 
derived by integration by parts [5] or other types of identities. One of the themes 
of today's computational technologies for evaluating master integrals concerns how to 
reduce a class of transcendental functions of the fornQ 

F(ci, . . . , cp; dx, . . . , d Q ; x x , . . . , x N ) 

ST^ r(ai-re + ci)r(a 2 -n + c 2 ) • • • T(a P -n + c P ) m nN 

2^ r(b v n + d 1 )T(b 2 -n + d 2 )---T(b Q -n + d Q ) Xl "' Xn U 

with dj-n = a^\n\ + • — \-a>i,N n N, etc. Generally, Cj's and d^s depend linearly on a small 
expansion parameter e, and the expansion coefficients in e need to be computed. The 
goal of the reduction is an expression in terms of some simpler mathematical objects 
whose properties are well known and which are amenable to fast high-precision numer- 
ical evaluations. Mathematical functions that have been explored for this purpose are 
the harmonic sum, nested harmonic sum |9J, polylogarithm, harmonic polylogarithms 
(HPLs) [10], multiple polylogarithm [IT], and their generalizations [12], [13] ; special val- 
ues of these functions, such as multiple zeta values (MZVs) [14] and their generalizations 
[T5| [12], [13] , also take part. 

The problem becomes more and more intricate as the order of loops is raised or as 
the number of scales involved increases. Typically in eq. ([TJ), the summation multiplicity 
iV becomes large, more T functions appear, more entangled combinations of the indices 
rii, ... , n at appear in the arguments of T functions, and expansions of q's and d^s around 
more complicated rational numbers and to higher orders in e become requisite. 

There have been systematic studies on reduction of single sums (N = 1) of the form of 
eq. ([!]) [T6] fT2 l [T7 1 [T5 | [I^ l l2T)]. By contrast, studies on reduction of multiple sums (N > 2) 
are still premature: The methods of recursion relations [21], of differential equations [22J, 
and of summation [T2T I23] have been examined and applied to computations of multi-loop 
integrals. What class of multiple sums can be reduced? And to which extent can they 
be reduced? Which method is more efficient or general? — Studies on these questions 
seem to have a long way ahead. 

In this paper, we present algorithms to evaluate two types of multiple sums, each of 
which corresponds to an arbitrary expansion coefficient in e of a sub-class of the above 
multiple sum. Type I sum corresponds to the case where the ratio of the T functions 
in the summand of eq. ([TJ reduces to a rational function of n^'s in the limit e — > 0, 

* For instance, this type of sums appear in the method of Mellin-Barnes integral representations for 
Feynman parameter integrals, after closing integral contours [2]. 
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with integer coefficients of n/s in each linear polynomial factor Furthermore, Cj, dj's 
can depend on an external integer index v linearly, through which the sum becomes 
dependent on v. The algorithm can reduce the expansion coefficients at any order in 
e and for any summation multiplicity. The result of reduction is expressed in terms 
of Z-sums (generalized nested harmonic sums) and their values at infinity (multiple 
polylogarithms), which depend on xi, . . . , x^. 

Type II sum is a double sum with the form £~ m=0 Up U q U r IfilZTuV 

where a p , . . . , f r are half-integers (including integers), and for p,q,r > 2, the ratios of T 
functions reduce to rational functions of £, m as e — > 0. This means that at most six T 
functions remain as e — > 0. The algorithm can reduce the expansion coefficients at any 
order in e. The result of reduction is expressed in terms of generalized MZVs. We give 
precise definitions for the two types of sums in the main body of the paper. 

Each type of sum covers a broad class of multiple sums and has useful applications. 
In particular, we can use Type I sum in the following ways: (a) to convert a non-nested 
integral to a combination of nested integrals, and (b) to derive a new type of relations 
among the generalized MZVs, which can be useful in reducing expressions of final results. 
We also apply the algorithms to compute master integrals for a static QCD potential at 
three-loop order. We present two new results. 

Our algorithm for evaluating Type I sum is similar to the one developed in [22] 
(see also [21]), which is an algorithm adapted to evaluating a wide class of multiple 
sums for loop computations. The algorithm of [23] uses solutions to difference equations 
for transforming a multiple sum to a combination of nested sums (which may later be 
expressed by harmonic sums, multiple polylogarithms, etc.) The method works if there 
are solutions to difference equations in the solution space in which they are searched 
for. At present, however, in many cases one has to make explicit trials to see which 
class of sums can eventually be transformed to nested sums. Specifically, it has been 
unknown whether Type I sum can be converted to nested sums in general cases. On the 
other hand, in our algorithm, we convert Type I sum to nested sums using the method 
of differences combined with the invariance of the summand under shifts of indices, 
in a recursive manner. This corresponds to constructing explicitly the solution to the 
difference equations with appropriate initial conditions, in terms of Z-sums and multiple 
polylogarithms, in the general case. Thus, our algorithm provides an existence proof of 
the solutions to the difference equations for Type I sums in general. Furthermore, our 
algorithm evaluates the sums fairly efficiently. 

Type II sums can be regarded as generalizations of the sums considered in [221 I2E1 
[T71 [20] . for the special values corresponding to Xj — 1. These previous studies evaluate 
expansion coefficients of similar transcendental sums around half-integer values of their 
arguments (not necessarily for Xj = 1), but they consider the cases where at most one 
pair of T functions remains in the summand after expansion (binomial sums or inverse- 
binomial sums). In contrast, we consider the cases where at most three pairs of T 



t This means that P = Q, all (cj — c?i) e _yo are integers, and ajj, bij as well as Cj, dj| e _>o are rational 
numbers. 
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functions remain, whose arguments include entanglement of the indices essential to a 
double sum. This adds significant complexity for evaluation of the sums in comparison 
to the previous works. 

In higher-order loop computations, more and more complicated classes of transcen- 
dental functions (and their special values) appear. In view of this general tendency, we 
use somewhat loose terminology in referring to generalizations of MZVs and HPLs, such 
as colored MZVs, cyclotomic numbers, etc. We simply refer to them as generalized MZVs 
or simply MZVs hereafter, and in the case that we focus particularly on their functional 
dependence on f3 1 (the upper bound of the outermost integral in integral representation), 
we refer to them as (generalized) HPLs, although they are essentially the same quanti- 
ties; they also belong to a certain class of multiple polylogarithms. See App. A for the 
definitions. 

The paper is organized as follows. In Sec. 2, we present the algorithm for reduction 
of Type I sums (Algorithm I). Some applications of this algorithm are shown in Sec. 3. 
We present the algorithm for reduction of Type II sums (Algorithm II) in Sec. 4. Sec. 5 
provides applications of Algorithm II. Summary and discussion are given in Sec. 6. We 
deliver definitions and conventions in App. A. Some details of Algorithm II are explained 
in App. B. 



2 Algorithm I: Multiple Sum without T Functions 

In this section we present an algorithm (Algorithm I) to reduce an iV-ple sum of the 
following formEl 

ii=a\{y) l2=a2(f,li) in=cin(l>,ii,...,in-i) 

where v is an integer. This includes, as a special case, in which S does not depend on v. 
In eq. @, A„ 6 C; each factor in the denominator L r (i>, i\, . . . , ijy) is a linear polynomial 
of v and i 1; . . . ,zjy with integer coefficients^; their powers p r E N satisfy ^ r p r > N; 
the upper or lower bound in each summation, a n or b n , is either infinity or a linear 
polynomial of its arguments with integer coefficients. We assume that, if v is within an 
appropriate range, the sum Siv) is finite, namely, the summand is finite for any values 
of the indices and the sum converges if some of a n , b n are infinity. After the reduction 
S(u) is expressed in terms of simple nested sums [generalized MZVs and Z-sums]. 



* Note that polygamma functions, originating from expansions of T functions in eq. ([T]), can be 
expressed using infinite series of rational functions. 

t Throughout this paper, "a linear polynomial of x\, . . . , x n with integer coefficients" represents 
cq + c±xi + ■ ■ ■ + c n x n with Cj € Z. 
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For clarity, we give simple examples]*] 

1 



oo oo 



7-r 7777 — ~ ; r = 4-2Cvr- — -41og2 + — C(3), 



oo k 



"I) 



k+m 



S£(* + 1) 2 (2™ + 1) 



-8Im[Z(oo;2,l;-i,i)] + 4Clog2 



7T 3 7T 2 

~8 ~ 12' 



(3) 
(4) 



n=0 m=0 v 



1 



n + is) 2 (l + m + n + z/)(m + n + 1) 



3; 1) + -Z{v- 1, 2; 1, 1) - ^Z(z/; 1; 1) 



(5) 



where C = YlT=o (2fc+i) 2 denotes the Catalan constant; ((z) = Y^Li 7^ denotes the Rie- 
mann zeta function; Z(oo; a%,--- , a N ; 6 1( • • • , an d -^(i/; a l, 4 4 4 , %; • • • , &jv) rep- 
resent a (generalized) MZV and Z-sum, respectively (see App. [X]for definitions). 
It is convenient to use the generalized sum defined by[§ 



£'/(*) 



k=a 



(o<6) 



o-l 



(a = 6 + 1) 
£ /(*) (a>6 + 2) 



fc=b+l 



for integers a and 6. The relations 

a 

£'/(*) = -£'/(*) 



6-1 



k=b 



k=b+/3 



k=a+l 
b+/3 



a+a— 1 



£'/(*) = £'/(*) + £'/(*)- £'/(*) 



fc=a+a 



k=a 



k=b+l 



k=a 



hold for arbitrary integers a, b, a and (3. Throughout this paper, it is understood that 
stands for . 



* The second example can be expressed also in terms of polylogarithm Li n (x) = J^feLi %a y i & 

7T 3 1 



\m[Z(pc; 2, 1; —i, i)] = — Im 



Li 3 



128 32 



7T log z 2. 



§The generalization is convenient for relating sums to functions such as harmonic sums and 
polygamma functions. For example, ^2 k=o k+T/2 = ^i+n l°g4 holds for an arbitrary integer n. 
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Let us define an integer matrix A from the coefficients of L r (u, i\, . . . , ijv)'s by 



Without loss of generality, we may assume that the rank of A equals the number of 
the summation indices N . In fact, if rank(A) < N, there exist redundant summation 
indices, namely, with an appropriate redefinition of summation indices, the summand can 
be rendered independent of N — rank(yl) indices, and the summation over these indices 
can be taken trivially. 

Algorithm I consists of the following steps: 

1. Introduce regularization parameters if necessary and perform partial fraction de- 
composition of the summand. 

2. Convert a multiple sum to a combination of nested sums. This is achieved by the 
method of differences for taking a sum of series, combined with appropriate shifts 
of the summation indices. 

3. Remove regularization parameters and convert the nested sums to MZVs and Z- 
sums. 

4. Reduce MZVs and Z-sums to simpler ones. 
We explain each step. 



Step 1 



We apply partial fraction decomposition to the summand of S(u) with respect to each 
summation index i n , starting from the innermost sum, up to the outermost sum. By this 
procedure, S(i>) is rendered to a combination of the sums of the form eq. ([2]) (apart from 
coefficients independent of i ly . . . ,i N ), with N factors in the denominator, L±, . . . , L N . 
Here, Li is independent of i 2 , ■ ■ ■ , in, L 2 is independent of i 3 , . . . , ijv, . . . , and L N _i 
is independent of i^. The matrix A is therefore rendered to an iV x iV lower triangular 
matrix. 

The partial fraction decompositions may generate divergent terms. By way of exam- 
ple: 



E 



1 



(*1 + *a)(*l +^3)(«2 +«3) 



V -1 

^ 2ii 



+ 



n + 12 n - 12 



i>2 + «3 H + H 



■ (7) 



The term — i 2 ) is divergent when %\ = i 2 . In general cases, if divergent terms are 
generated, we introduce regularization parameters 5 n and shift0 i n — > i n + 5 n . At a later 



™ The advantage of this regularization is that the limits i\ — > «2 and Si, 62 — > commute. 
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stage of the computation (step 3), we expand the result of the computation in 5 n . All 
singular terms in 5 n will cancel outH 

Step 2 

After step 1, we obtain a combination of sums, each of which has a form 



s,m= e e ••• e n . r \ < 8 > 



U=ai(^) ii=a,2{vj,\) i N =a N (u,h,...,i N _ 1 ) 



up to an overall coefficient independent of ii, . . . ,i N . Here, (5 = (Si, ... , 5n). The goal 
of step 2 is to convert S 2 to a combination of nested sums; a nested sum has a general 
form 

E Atii,8)f2(j2,S)---f M (3M,S), (9) 

A(^)>ii>->iM>i 

where A(i/) is either oo or a linear polynomial of v with integer coefficients. We construct 
an algorithm by induction. Suppose that for multiplicity iV < K — 1 we have a procedure 
to convert the sum given by eq. (jSJ) to a combination of nested sums. Below we show 
how to convert the sum with multiplicity N = K. 

We first find a set of integers (Ao, Ai, . . . , An) which satisfy the following two con- 
ditions: 

(i) For all L r , 

L r (u + A , z'x + Ai, . . . , i N + A N ; S) = L r (v, i u . . . , i N ; 6). (10) 

(ii) A ^ 0. 

From the condition (i), it follows that ^ n ^4 rn A n + C r A = 0, where C r = dL r jdv. 
Hence, 

A n =-J2(A- 1 )nrC r A . (11) 

r 

Since (A~ l ) nr are rational numbers, with an appropriate integer A 7^ 0, all A n can be 
set to integers. We denote 

(v + A ,Zl + Ai, . . .,i n -l + A n _i) - an(v,ii, ■ ■ -,in-l), (12) 

Pn = b n (v + A , ii + Ai, . . . , i n _! + A n _i) - b n (u, h,..., i n -i), (13) 
which are also integers. 



II For practical implementation, it is often more efficient to introduce a single regularization parameter 
5 and to shift i n y i n + c n S with an appropriate choice of constants c n s. 
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Using the property (i), one can relate S 2 {y + A ) to S 2 (v). In fact, if we define the 
"difference" of the two terms as 

AS 2 (u) = S 2 (u + A ) - \S 2 (u) ; \ = l[\^, (14) 

n 

we may rewrite 



AS 2 M = A 



bi+f3i b N +f) N bx b N 

E -. E -E-E 

ii=ai+ai ijv =a iV+ Q iV ii=<n *jv =a iv. 



[] r L r (is,ii, ...,i N ;S) 



Pr 



(15) 



by shifting the indices as z„ — )■ i n + A n in ^(z/ + A ). With decomposition of the sumo 

bi+fii b N +/3 N bx b N 

E • E -E-E 

61 ai+ai— 1 \ / b N a N +a N -l b N -\-{3 N \ ^ 

E- E E ■■■ E - E + E -E - E. 

(16) 

bulk of the difference cancels. The remaining terms include at least one "surface term," 

sum can be evaluated explicitly since a n , f3 n are explicit 
numbers. Thus, the remaining terms have summation multiplicity K — 1 or less. Ac- 
cording to the assumption of induction, AS 2 (u) can be expressed by a combination of 
nested sums, eq. (jUJ). 

We can reconstruct S 2 (v) by summing up an initial term S 2 {vq) and the differences 
AS 2 (v) with appropriate weights. Here, vq is the remainder in division of v by Aq, i.e., 
v = vq (mod Aq) with < uq < A . We have 

u-A 

S 2 {p) = \^ )/A °S 2 (u ) + Proj(^,^o, A ) ASM A^- A °)/ A °. (17) 

The projector is defined as 

„ , 1 / , m — n\ \ if m ^ n (mod d) 

Proj m, n, d = - > exp ( 2nik — ) = < ^ , 18 

JV ; d £i P V d J I 1 if m = n (mod d) V ; 



fl=U 



k=l 



such that the sum in eq. (TTTj) is taken in steps A . Due to the presence of the projector, 
the upper and lower limits of the summation in the second term of eq. ( ITT]) can be 
replaced by v — 1 and 0, respectively. 



** The decomposition may generate divergent terms. In this case, one needs to introduce regularization 
parameters as in step 1 (if not introduced already). 
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In the first term of eq. (fTTl) . if we assign a specific value to z/ G {0, 1, 2, . . . , A — 1}, 
^(^o) can be converted to a nested sum. In fact, denoting ^(z/q) = ^ -F^i), we may 
regard %\ as an external index of a sum -F(h) of multiplicity K — 1; hence, -F(^i) can be 
converted to a combination of nested sums. 

Since u is a function of v in eq. (JTTJ) , for our purpose, we rewrite the right-hand side 
in the following manner: we replace vq by /io, multiply by a projector which projects out 
the term fio = uq, and take a sum over /io- Hence, 



An-l 



S 2 ( U ) = \^° ^Proj(z/,/i ,A c 



x 



M>=0 



A-^ A "S 2 (/i ) + ]T Proj(/i, /i , A ) A5 2 ( M ) A-^ +A °)/ Ao 

/u=0 



(19) 



A being an explicit number, the sum over /i can be evaluated explicitly, and the result 
is A V//A ° times a combination of nested sums. Thus, apart from v dependent coefficients 
which multiply individual terms, S^iy) is converted to a combination of nested sums. 
Let us present a simple example: 

00 {—\\ k 

f( m ) = jT^rrh^- ( 2 °) 

' (1 + k + mr 

fe=i ' 

This can be expressed as 

f(m) = [f(m) + /(m - 1)] - [/(m - 1) + f(m - 2)] 
+ - + (-l) m [/(2)+/(l)]-(-l) m /(l) 

m 

= E(- 1 ) m_i [/0') + /0' - !)1 - ("l) m /(l). (21) 

3=2 

Bulk of the sum in the "difference" of two adjacent terms gets canceled since the shifts 
m — > m — 1, k — > k + 1 leave the denominator in eq. (120]) unchanged: 

/ oo oo \ / -.xfc _ 1 

/0) + /o-i) = (E-E)(y F =(-^ T? - (*> 



Thus, 



/(») = (-i>-E^-(-i>-E^- < 23 > 

In the first term, apart from the coefficient (— l) m , dependence on the external index m 
enters only through the upper bound of the summation, while the second term, apart 
from (— l) m , is free of m and is essentially an MZV. 
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The above relation can be used to convert a double sum without an external index 
to a combination of nested sums: 



y^y> (-l) k i m im /( m ) 

^ ^ (2 + m) 2 (l + k + m) 2 ~ ^ (2 + m) 2 

(2 + m) 2 ^ (j + If 2^ i i2 + m f^ [ {k + 2r 1 J 

The first term in the last line is essentially a nested sum, while the second term is a 
product of single sums. 

Step 3 

We convert the nested sums obtained in step 2 to MZVs and Z-sums. Up to an overall 
coefficient independent of the summation indices, each nested sum has a form 



S 3 (u;S)= Y ^ X--.X- ^ : (21) ) 

A(f)>Ji>— >JM>0 

where c n = (c nj i, . . . , c^at) and c njr , g n G Q. We have normalized the coefficient of each 
summation index j n to unity in each factor in the denominator. 

First we convert all the offsets g ra 's to integers. Let £ be the least common multiple 
(LCM) of the denominators of qi, . . . , qu- In eq. ( )25|) . we may replace j n by j n /£ and 
insert Proj(j n , 0, £) for each n, such that each sum is taken on multiples of £: 

^A(^)>ii>->j A />0 n=l VJn yn n ; 

Expanding the projector eq. (fl8l) . S3 is given by a combination of nested sums with 
integer offsets £q n . 

Next we shift each index j n — > j n — £q n to eliminate the offset. By this procedure, 
divergent terms as 5 n — > are explicitly taken outside of the summation. Then we 
expand in 5 n 's [up to C(<5°)] the terms taken outside of the summation as well as the 
summands. All the divergent terms (which include negative powers of <5 n 's) cancel out. 
Finally, by repeated applications of shifts of summation indices, shifts of upper and 
lower bounds of sums, and partial fraction decompositions, we can convert all the sums 
to Z-sums [with argument nv (n G N)] and MZVs. 

In the case that some of A(i/)'s are 00, divergent MZVs as j n — > 00 may appear 
in individual sums. By regularizing the sums carefully for large values of j n 's, these 
divergences can be extracted from individual terms and they cancel out. For instance, 
if we adopt a cut-off regularization, A(z/) = 00 — > A(z/) = A 3> 1, at most logarithmic 
divergences ~ (logA) p appear. A shift of index j n alters neither divergent part or the 
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finite part as A — > oo. A multiplication of the cut-off in eq. (126]) induces a change 
(log A) p — > (log A + \og£) p , hence, this effect must properly be taken into account. 
We give an example to illustrate the step 3: 

y- 1 y- 1 = y^ Proj(fc,0,2) y- Proj(m, 0, 2) 

2s k — 1 -I- 5 2-^1 m .4-1 /9 L{./9-UaL 



fc=l 

2v 



k-l + 5 ^ m + l/2 ^k/2-1 + 5^ m/2 + 1/2 

m=l ' fe=l ' m=l ' ' 

l + + 2 ^l + (-lf Al + (-l) 



fc-2 + 25^ m + 1 35 ^ - 2 + 25 ^ m + 1 

fc=l m=l fc=3 m=l 

2 52 2 4 Z(2i/; 1; 1) - Z{2v- 1; -1) 



35 9 3v 3(1 + 2i/) i/ 

- 2Z(2i/; 1; 1) + y Z(2i/; 1; -1) + Z(2^; 1, 1; 1, 1) - Z(2u; 1, 1; 1, -1) 

+ Z(2z/; 1, 1; -1, 1) - Z(2u; 1, 1; -1, -1) + 0(5). (27) 

In the last equality, a straightforward but somewhat cumbersome computation is needed 
to convert the sum to a combination of Z-sums. 



Step 4 

In many cases, host of MZVs and Z-sums are included in the result of the step 3. They 
can be expressed by small sets of basis elements for MZVs and Z-sums through specific 
reduction procedures or by utilizing known database. Reductions of MZVs and Z-sums 
using various relations among them, such as shuffle relations, have been studied exten- 
sively in the literatures, and we do not discuss them here; see [271 [13] and references 
therein. We apply such reduction procedures to simplify the result. (We can also derive 
non-trivial relations among MZVs using the algorithm given in this section. This will be 
discussed in Sec. 13.21 ) 



Comment on surface terms at infinity 

In the case that some of the upper bounds of the summation in eq. (j2J) are infinity, the 
partial fraction decompositions in step 1 may generate logarithmically divergent sums 
as jn —> oo. This causes some subtleties to contributions of the "surface terms," which 
appear in step 2. For convergent sums, the surface terms from j n — > oo certainly vanish. 
It is a priori not obvious, however, whether contributions of the surface terms vanish 
in the case that there are individually divergent sums. By way of example, consider a 
convergent sum (=1): 

k(k + n)(k + n + 1) ~ ^ ^ k \k + n k + n + 1 J ' ^ 

k=l n=l \ / \ / k=l n=l x ' 
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If we separate the sum on the right-hand side, each sum becomes logarithmically diver- 
gent. Let us replace the upper bounds by a cut-off A ^> 1 and convert the first sum 

ELi with f( k ) = T,i=i fci^ to a nested sum b y 



/(*) = /(I) + £[/(?) " M - !)1 = /(!) + E ("TA ~ ~) 
i=2 i=2 \ J + J / 

Thus, the contribution of the surface term at j — A is given by 

r 2 



(29) 



El 1 7T 

fc ^ j + A ~ 12 

=1 3=2 J 



as A oo, (30) 



k=i j-~ 

which is indeed non-vanishing. It turns out that the contribution of the surface term 
from the second term of eq. (|28|) exactly cancels eq. (|30|) . Hence, the contributions of the 
surface terms vanish as a whole. This would not be surprising, since the sum eq. ( 128|) is 
convergent as k — > oo. Instead, if we introduce different cut-offs A^ and A n for k and n, 
respectively, and take the limit A n — > oo first, we readily find that the individual surface 
terms vanish. 

In general, this problem of the surface terms can be treated properly by introducing 
cut-off regularization as above. In a particular regularization scheme, we can argue that 
the surface terms from j n = oo vanish (provided the computation is carried out in a 
specific order, see below), so that these surface terms can be ignored in the computation. 
This reduces the amount of computation considerably, as compared to other regulariza- 
tion schemes, in which one has to trace contributions of the surface terms through the 
computation. 

The argument goes as follows. Suppose that, in the original sum S [eq. ©j, the 
upper bounds of the summation indices i nr for r = 1,2, . . . , N'(< N) are infinity. We 
replace the upper bounds of these indices by A Hr , respectively. S is convergent as each 
A„ r — > oo. For convenience, let us denote the external index v as iq. The surface terms 
are contained in AS2, eq. ffl5l) . Since we apply the algorithm in step 2 recursively to 
convert sums with lower multiplicities, the argument i x of 5*2 or AS2 is not necessarily 
an external index (i ) but generally it can also be one of the summation indices i n 's. 
We can always work from the outermost to innermost sum in applying the algorithm, 
such that any of the indices corresponding to the surface terms in AS2, i n , is inner with 
respect to the argument i x of AS2, namely, x < n for all the summation indices z„'s 
contained in S2 (i x )- Consider the contribution of only one surface term at i Ur = A„ r to 
AS 2 : it is suppressed by with p > 1; contributions of sums of indices other than i Hr 
diverge at most logarithmically ~ Yl s ^r(^°S -^■n 3 ) Ps ■ If the index i x reaches up to order 
A x , the sum of AS , 2 (2 X ) over i x is at most order A^/A^ x Y[ s ^r^°S-^ns) Ps ■ Therefore, 
if we parametrize all the cut-offs by a single parameter u as A„ r = u nr and take a limit 
u — > 00, the cut-offs corresponding to inner indices diverge faster (A x A„ r ), so that the 
contribution of the surface term vanishes. (The outer sums modify the dependence on u 
at most logarithmically, which does not alter vanishing of the surface term contribution.) 
The contribution of more than one surface term vanishes even faster with u. 
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To summarize, if we work from the outermost sum to innermost sum in converting 
non-nested sums to combinations of nested sums (if we do not exchange the order of 
summation indices), we can neglect contributions of the surface terms. 



3 Applications of Algorithm I 

3.1 Conversion of Non-nested Integral to Nested Integrals 

The algorithm given in the previous section (Algorithm I) converts a certain type of 
non-nested sums to combinations of nested sums. A straightforward application of this 
algorithm is to convert a certain type of non-nested integrals to a combination of nested 
integrals such as HPLs. In fact, in problems suited to manipulations in integral repre- 
sentations, it is often useful to convert a non-nested multiple integral to a combination 
of nested integrals. 

There are multiple integrals (with one external parameter x) which can be converted 
to multiple sums of the form 



F(x) = J2^S(v)Proj(v,0,d) ; PeN, 



(31) 



by first expanding the integrand in Taylor series and then integrating. Here, S(u) is 
given by eq. fl2]). This may occur in the case that the integrand is a combination of a 
rational function, logarithms, and HPLs of the integral variables. We can use Algorithm 
I to convert S(u) to a combination of nested sums; in this form F(x) is expressed 
as a combination of nested integrals straightforwardly. We give a simple example for 
illustration: 



Kx) 



dy log(l - y/x) log(l - xy) 



k=l 



(xy) 1 



m 



00 x l+2m 00 

m 



x x l+2m 



m=l 



-^k(k + m + l) 



£ 

m=l 



m 



m=l 

Z(m;l;l) 



(m+1) 2 



+ 



m + 1 



(32) 



In the last step we applied Algorithm I. Using standard relations among MZVs and 
HPLs, we may express the above result in terms of HPLs: 



Kx) 



2x - 
1 



1-x 2 



2x + 



x 



X 



Z(oo; 1; x 2 ) + Z(oo; 2; x 2 ) + Z(oo; 1, 1; x 2 , 1) 



x 



H-i(x) + iZi(a;) - - H- 1A (x) 

+ 2iJ ,-i(x) + 2H 0>1 (x) - - Hi tl {x) 



(33) 



This is a combination of nested integrals; see eqs. ()671)-( l69i) for the definition of HPL. 
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In more general cases, one may need to redefine the summation indices appropriately, 
so that the exponent of x is expressed in terms of a single index. For instance, we redefine 
k + n = I in 



oo I— 1 



k=l n=l 1=2 n=l 



and convert the inner sum to a combination of Z-sums with arguments mi (m G N). 

Examples of integrals, which can be converted to combinations of nested integrals in 
a similar manner, and which may be useful, are given by 



h(x) = / dyy^ (\ ogy yW [log(1 _ y)] p(3) [log(1 + y) ]P(4) 

x [log(l - y/x)]*® [log(l - zy)]^ 6 ) tf(y), (35) 

7 2 H= / dx dy dzx p Wy p( ®z p W [log (1 -xy/w)] m 
Jo Jo Jo 

x [log (1 - yz/w)f 5) [log (1 - zx/w)] m , (36) 

where p(n) G N, and H(y) is an HPL. The resulting expressions, however, may be quite 
lengthy. 

3.2 Relations among MZVs 

As already stated, a number of relations among MZVs have been derived and used for 
expressing MZVs by a set of basis elements. We show that Algorithm I can be used to 
derive yet another type of relations among (generalized) MZVs. 

We demonstrate it using an example. Suppose that A G C satisfies 

|A| < 1 and A ^ 1. (37) 

We rearrange the order of the summation of a weight-two MZV: 

L z J ( -\\m. 



Z(oo;l,l;A,-A)= £ = £ £ A - 

nm ' mis — m) 

n>m>0 s=3 rn=l 

oo r , , m oo r / 1 \m 

= v \ 2r+i y — ^ — - + v A 2r+2 v /o ( ~^ — -. (38) 

^ m(2r + l-m) ^ ^ m(2r + 2 - m) v ' 

r=l m=l v ' r=l m=l v ' 

We set s = n+m. The upper bound of m follows from the condition m < n—\ = s—m—1, 
where |_^J denotes the greatest integer less than or equal to x. In the last line we separate 
s to odd and even numbers and set s = 2r + 1 and s = 2r + 2, respectively. The sums 
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over m can be converted to combinations of Z-sums by Algorithm I. It can be shown 
that the above rearrangement of summation is justified under the condition eq. fl37|) . 
Thus, we find a relation among MZVs: 



z(X, -A) = --z(-2X 2 ) + -z(-X, -1) - z(-X, -i) - z(-X, i) - -z(-X, 1) 

- \z{X, -1) + z(A, -i) + z(A, i) - \z{X, 1) + - A z{X\ 1). (39) 

We use a short-hand notation for MZV, eq. ( 166|) . where the original forms can be re- 
produced via z(x\, . . . , x n ) Z(oo; a 1; . . . , a n ; 0i, . . . , (3 n ) with = |xj(A = 1)| and 
A = Xi/aci- The relation eq. ( ]3"9~]) may be regarded as a special case of the conversion 
relation described in Sec. 13.11 

In particular, in the case that A ra = ±z for some n e N, eq. ( 139]) may lead to an 
interesting relation. We have examined independence of the above relation from other 
relations, such as shuffle relations and those which can be derived by variable transfor- 
mations in integral representations of MZVs (e.g. x — > x , x — > x 2 , x — > 
the case that A = ±z, eq. (13"9"|) does not give a new relation. On the other hand, in the 
case A = e™^ (a primitive 8th root of unity), it gives a new relation, which can be used 
for reduction of MZVs. In fact, this gives a relation between the three basis elements at 
w = 2, / = 8 proposed in [T3] : 

5371-2 

-Si Im[Li 2 (e 8 )] + 3a hl (e 8 , e 4 ) + a hl (e 8 , e 3 8 ) + 2iC + 

+ ilog 2 (72-l) - ^ - \ log 2 log^-l) + ^vrlog2 = 0. (40) 

Therefore, the number of the (new) basis elements reduces to two. 

It is easy to obtain relations among MZVs in more complicated cases, by rear- 
ranging summation orders using this method. The strategy is to rewrite an MZV 
Z(oo; ai, . . . , a n ; (3i, . . . , (3 n ), for which some of are the same (equal to A), in a com- 
bination of HPLs of the type i? aiv .. )an (a;A m ) (MZVs with (3\ = uX m ). Here, u denotes a 
root of unity. 

Up to now, we do not know whether these relations can be derived from other known 
methods. Extensive exploration of the relations among MZVs for the case beyond Pi = 
±1 (e.g. n-th root of unity for n > 2) is still underway [151 1211 [13]. We expect that the 
above method would be a useful tool for analyses in this direction. We have checked 



* Using these relations we obtain, for instance, a relation between the basis elements at w = 3 for 
the infinite cyclotomic harmonic sums over {(±l) fc /fc, (±l) fc /(2fc + 1)} proposed in [13] : 

1 7T 3 7T 2 

-JO"{1,0-2},{2,1 -1} + (7 {2,i ,-2},{i,o,-i} +Clog2 - — + — - 1 + log2 = 0, 

hence, ((3) and Im[Li3(i±i)] (= —C log 2 — ^ct{ 2 , 1,-2}, {1,0,-1} + ffg- + §> l°g 2 2) suffice to be introduced 
at this weight. 
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k + q 



k + q" 



k + K 




P 



Figure 1: Diagram for the master integral Double lines represent heavy-quark propaga- 
tors, while single lines represent massless propagators. 

in specific examples that they lead to relations independent of shuffle relations and are 
powerful in reducing MZVs to a small set of elements. Generally, a huge set of relations 
need to be generated, which requires an efficient algorithm for generation of the relations, 
and Algorithm I is useful in this respect. 

3.3 Evaluation of a 3-Loop Master Integral 

We apply Algorithm I for evaluating a 3-loop integral 



This is one of the 40 master integrals, which appear in a computation of the 3-loop 
correction to the static QCD potential (a 3 ) [221 ED]; the diagram is shown in Fig. [TJ In 
dimensional regularization (D = 4 — 2e), expansion coefficients of Ai\ in e up to (and 
including) order e is necessary. In eq. ( 14T|) . we omit +i0 in the propagator denominators, 
hence, {k-v , k 2 , (k + q) 2 , . . . } stand for {k-v + iO, k 2 + i0, (k + q) 2 + i0, . . . }; q — (0, q) is a 
spacelike external momentum, whereas v = (1, 0) is a temporal unit vector; {k ■ v + iO)^ 1 
represents a heavy-quark propagator. 

Using a standard technology [2], such as the one described in [31], one can convert the 
above integral to a three-fold Mellin-Barnes integral. Closing the contour in the complex 
plane for each Mellin-Barnes integral variable, one may further convert the integral to a 
combination of three-fold sums with T functions. By expanding in e, at every order of 
e all the T functions cancel, and we are left with sums of the form eq. ([2]). We evaluate 
the expansion coefficients up to 0(e) using Algorithm I, where sums with multiplicity 




(2%) D (2n) D (2rr) D (k ■ v)(p ■ v)k 2 p 2 (k + q) 2 (p + q) 2 (k + k) 2 ( P + «) 



(41) 
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up to six need to be evaluated. The result is given by 

^4 



M 1 = A(e)(-q 



2\-l-3e 



71 

e 



-93C(5) - 14tt 2 C(3) + tt 4 ^2 + 6 log 2) J 
+ e |-288s 6 - 186C(5) + 234C(3) 2 + 7r 2 ^96Li 4 (|) - 28C(3) + 4 log 4 2 

+ 7 r 4 (4+121og2 + 141og 2 2)+^} 



A(e) = -i (47r) 3e - 6 exp(-37 i; e 



(43) 



(-1)" 



where je = 0.5772... denotes the Euler constant; sq = z(—5, — 1)+£(6) = ^2 m>n>1 m5n — 

0.9874414264...; Li n (x) = J2T=i fn denotes polylogarithm. We have checked the above 
result numerically by evaluating Feynman parameter integral representation using sector 
decomposition. This result has not been published elsewhere. 



4 Algorithm II: Expansion Coefficients of Double 
Sum with r Functions 

We consider a double sum 

00 00 

K(a,b; c,d;e,f) = ^ 

e=o m=0 

(44) 

where a = (01, . . . , a^), etc. In this section we present an algorithm (Algorithm II) for 
computing expansion coefficients of K in a small parameter around integer or half-integer 
arguments. We assume the following conditions for the arguments of K: 

(i) a = a + eat, ■■■,/ = fo + e /i> where a , p , b 0>p , c , p , d , p , e 0tP , f 0>r are half-integers 
(including integers) and e is a small expansion parameter. a% !P , b% iP , ci tP , d\ tP , ei tP , fi >r 
are complex variables used to parametrize the expansion coefficients. 

(ii) For p,q,r> 2, a , p - b 0iP , c 0j(? - d 0>q , e 0;r - f 0jT are integers. 

(iii) The sum in eq. ( 1441) is convergent at (a , bo] Co, d ; eo, f ), hence K is regular with 
respect to 01, . . . , / x as e oE 

The algorithm reduces an arbitrary expansion coefficient of K in e, up to any order, to a 
combination of MZVs. Due to the above conditions, there remain at most six uncanceled 
T functions in the summand in each expansion coefficient, before reduction. 

* It is important that K does not include a singularity such as ^xg^ around e = 0. 
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nr(e + a p ) 
r 



M 

n 

L<7=1 



r (m + c q ) 
r (m + d a ) 



nr(£ + m + e r ) 
_ r(£ + m + / r ) 



We give an example: 



d d d 2 
de 85 da 2 



oc 



r(£; + l)r(m+l)r(£;+m+f)r(£;+m + e+2) 



C(e,8,a) 



r(A;+a + |)r(m+|)r(A;+m+2)r(A;+m + (5+3) 



k,m=0 



e,5,a— >0 




+ C(5) (6882 log 2 + 22318) + C(3) (-448 log 3 2 - 2024 log 2 2 - 752 log 2 - 248) 



1024 log 6 2 3808 log 5 2 32 log 4 2 1024 log 3 2 , , / nnnnT . , 1N 
jf jf f— + f— + log 2 2 (-3072Li 4 (I) + 192 



) 



+ log2^-9216Li 5 (i) -8320Li 4 (|) - 16o) 

- 12288Li 6 (|) - 11136Li 5 (|) - 256Li 4 (±) - 144, 



(45) 



where we have included a prefactor C(e,5,a) = r(<5+3)r(o; + |)/r(e+2) in order to 
eliminate trivial 7# and y/n. 

The algorithm consists of the following steps: 

1. Convert K to an integral representation. 

2. Apply Rummer's formula and other variable transformations to simplify the inte- 
gral and to eliminate half-integer powers in the integrand. 

3. Expand in e. If necessary, regularize appropriately before the expansion. 

4. Convert multiple integrals to combinations of nested integrals. Integrals with re- 
spect to individual integral variables are converted recursively. In each conversion 
process, we extract singularities using integration by parts, remove regularization 
parameters, and apply Algorithm I. In the end the result is expressed by MZVs. 

5. Reduce MZVs to simpler ones. 
We explain each step below. 

Step 1 

For later convenience, we denote 



e r = e' N+1 _ r , f r = f' N+1 _ r (1 < r < N). 



(46) 
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We express each pair of T functions using the integral representation of the beta function 
as 



r (l + dp) _ 1 f 1 J„ £+a p ~l (1 _ \b p -a p -l 

T(£ + b p ) T(b p -a p )J dXpX ? [i Xp) 



T(l + m + tt) T(fl-e' r )J 
Hence, we obtain an integral representation 



pi pi pi pi pi 

K oc / dx\--- / cLxl / dyi--- I dy M / dz\ • ■ • / dz N 
Jo Jo Jo Jo Jo 



(47) 



T(m + c q ) 1 r +Cq _ 1( _ , dq -c q -i (Ag) 

T(m + d q ) T(d q -c q )J ^ [l Vq) ' ^ 

1 (/ ' "' " 41 = 1 f dZr - Zr )fr-<-\ (49) 



X ' 



1 - aci • • • x L zi • • • - y x ■ ■ • y M z 1 ■ ■ ■ z N ) 

(50) 



Step 2 

We apply Kummer's formula 

r 1 t a {i-tf r(a + i)r(/3 + i) (i_t)^ +1 r T t a+ ^ +l ( t\~ l 

J 1 (1-tTp - r( 7 )r(a+/3- 7 + 2) T«+^+ 2 J * (1 - t)W \~Tj 

(51) 

for 7 = 1 to the variables t = X\ and y\ of eq. (150]) . Next we apply the following variable 
transformations successively: 

• zn = Cat, and Zi = Q/Ci+i starting from i = N — 1 to 2 = l[j] 

• Xi = £i, = £l/Ci, an d x « = 6/6+1 starting from i = L — 1 to 2 = 20 

• Ui = 2/m = and ^ = r/i/^+i starting from i = M- ltoi = 2. 

This leads to an integral representation of K, where the integral variables are ordered as 
depicted in Fig. |2j and where each factor of the integrand contains at most two (adjacent) 

t If JV = 1, Z! =Cl- 

* If L = 1, x\ = £i; if L = 2, x\ = £i and = ^/Ci- Similar transformations apply in the cases 
M = 1,2. 
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Figure 2: Diagram showing the order of the integral variables. Along the same line, the 
variables to the left are larger. The order of £ p and r\ q is not fixed. 



integral variables: 

K oc fd( N •••<!& [ #1 / d VM ■■■d Vl Cf^™- 1 (1 _ Cn/^- 1 
i>Civ>->Ci>o Ci>Ci>->5i>o Ci>»?m>->»?i>o 



x 1-- 



C \ b L -a L -l / \ d M ~c M -l 

Si 



p \ 62-012-I / t \ 6 £-l - a -L-l— 1 

S2 \ ^. a3 _ a2 _i / 1 SL— 1 \ t .a L -a L _ 1 — 1 



x(i-*) % — •••(1-^) 

SI \ >e'-e -1 / 1 C.JV-1 X 



x^i-gj c 2 • • • - c; ■ (52) 

We note that all the factors containing two variables have integer powers at e = 0, 
according to the condition (ii) and eq. (|46|) . On the other hand, the factors t and 1 — t 
for t = £ p , r] q , ( r have half-integer (including integer) powers at e = 0. 

We may render the exponents of all the factors of the integrand to integers (at e = 0) 
using the following variable transformations. We first apply the transformation 

t ->• — (53) 

to all the variables t = C, p ,r] q ,(r simultaneously. The integral region is unchanged, and 
the factors in the integrand transform as 

\l + t) V t'(l+i) 2 k ' 

If half-integer powers still remain, we further transform t — > t 2 . This eliminates half- 
integer powers in the general case. 
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Table 1: Exponents of the factors of the integrand in eq. (|52p after introducing regularization 
parameters eq. (f55jh The exponents of £1, £2, f?2, respectively, depend on two regularization 
parameters. All other exponents, respectively, depend on a single regularization parameter. 



Step 3 

Although K is finite at e = 0, setting e = in the integrand before integration may lead 
to divergences. These superficial divergences originate from the integral representations 
eqs. (|47l)-(l49l) and ( JoTT) . in the case that some of these integrals are not well-defined 
around e = 0. If there are (superficial) divergences, we introduce regularization param- 
eters in the following manner, before expanding in e: 

Oi -> a, + 8 ai H h 8 ai , h -4 bi + 8 ai H h 8 ai + 8 bi , 

Ci -)> Q + 5 Ca H h 5 Cl , di ^ di + 8 Cl H h 5 Ci + 5^, 

e- ->■ e- + 5 ai + • • • + 5 aL + 8 C1 + ■ ■ ■ + 8 CM + 8^ + ■ ■ ■ + 8 e >., 

fi-> f , i+8 ai + --- + Sa L +6 cl + --- + 5 CM + 8 e[ + ■ ■ ■ + 5 e > + 5 n . (55) 

This parametrization is specially tailored for the use in step 4. As we will see, it is 
important that we can take the limit 8i — > corresponding to each integral variable (each 
set of integral variables), without affecting outer integrals. In fact, after applying the 
above regularization, most exponents of the factors of the integrand depend on only one 
regularization parameter, respectively. For reference, in Tabled] we list the exponents in 
eq. fl52l) after applying the regularization, from which one may easily deduce the integral 
form after the variable transformation eq. fl53|) (and t — > t 2 ). 

We expand K in e (if necessary, after applying the above regularization). The expan- 
sion generates powers of logarithms [logt, log(l±£), log(l±it), log(l±i/t'), log(l±tt')] 
in the integrand. 



Step 4 

We convert the integral to a sum of nested integrals using a recursive algorithm. The 
order of conversion is indicated by the arrows in Fig. [2j we work from the innermost 
integral up to the outermost integral. First we convert the double integrals with respect 
to the innermost variables £1,^2 and 771,772, respectively. For the other variables, we 
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convert the integral with respect to each variable recursively!! 

Let us explain how to convert the double integral with respect to £1 and £2- The 
integrand depends on the four regularization parameters 5 ai ,5 a2 ,5b 1 ,5b 2 (denoted by 6^2). 
We note that these parameters are not included in the outer integrals. Using identities 
and integration by parts, we raise or decrease the exponents of the factors, £j, l±£j, l±i£j 
(i = 1,2), l±£ 2 /£3, 1±^2^3, until all of them become non-negative when ^12's are set 
to zero. This procedure is explained in App. [Bl By this procedure, the double integral 
is converted to a sum of double integrals, where singularities in ^12 's are taken outside 
of each double integral as poles and each double integral is finite as all 6^2 — > 0. We 
expand K in ^i 2 's up to (and including) finite terms. The poles in ^12's cancel, and 
each double integral with respect to £1 and £2 is reduced to the form 

#2 J d& P €12 (6,l0g6,l0g(l±^),l0g(l±^) ) l0g(l±6/e3),l0g(l±6£3)) 

(56) 

where P^u is a polynomial of its arguments. Prefactors dependent on £3 may be gen- 
erated, which are absorbed into the outer integral J 0^3. The above integral can be 
converted to a combination of HPLs using Algorithm I, as described in Sec. I3.1|P1 In 
fact, it is a generalized version of eq. (155]) . 

In the same way, we convert the double integral with respect to rji and r] 2 to a sum 
of HPLs. 

Next we apply the following algorithm recursively from inner to outer integrals. Sup- 
pose we convert the integral with respect to t(— £ p (>3), %(>3)> CO- As a result of the 
previous conversion, each integral has the following form: 

dt n^(MT (5l ' &) P(log0 1 (t,t , ),---,log0 9 (t,O)^W- (57) 

_s=l J 

Here, <fi s {t,t') denotes one of the nine factors t, (1 ± £), (1 ± it), (1 ± t/t'), (1 ± tt'). 
P(xi, . . . , xg) represents a polynomial of xi^...,x s . Hit) is a nested integral which 
resulted from processing the inner integrals!!!! In fact, Hit) is a generalized HPL of 
eqs. (!67|) - (l70|) . There are only two regularization parameters associated with the variable 
t, and we denote them by 61,62- In particular, 61,62 are not included in the outer 



§ The cases with L = 1 or M — 1 are exceptional. If L = 1 and M > 2, we work in the following 
order: we first convert the double integral with respect to 771 , 772 ; we work along the rj chain for each 
variable recursively up to r\M\ we convert the double integral with respect to £1 and £1; we work along 
the £ chain recursively up to £jv- (If M — 1 and L > 2, exchange and fy.) If L = M = 1 we convert 
the triple integral with respect to £1, rji, £1 first; the method is similar to the double integral case. 

^ All the logs in Pji2 can be expressed by an HPL. It is useful to use the shuffle relations (fusion 
rule) of HPLs to express their products as a sum of HPLs. 

" For t — (1, H(t) is replaced by a product of H^(t) and H v {t), which originate from £ and 77 chains, 
respectively; see Fig. [5] One can use the shuffle relations (fusion rule) to express H^it)H„it) as a sum 
of nested integrals. 
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Figure 3: Diagram for the master integral M.2- Notations are the same as in Fig. [TJ 

integrals. The exponent a s (5i,5 2 ) is a linear function of 61,62, where a s (0,0) is an 
integer. 

We raise or decrease the exponents a s (5i,5 2 )'s using identities and integration by 
parts (see App. [B|. By this procedure, the original integral is converted to a combination 
of integrals, where singularities in 61,62 are taken outside of each integral as poles and 
each integral is finite as 61,62 — > 0. We expand K in 61,62 up to (and including) finite 
terms. The poles in 61,62 cancel, and each integral with respect to t is reduced to the 
form 

J dtP(t,\ogMt,t'),...,\ogMt,t'))H(t), (58) 

where P is a polynomial of its arguments. Prefactors dependent on t' may be generated, 
which should be absorbed into the outer integral0 The above integral can be converted 
to a combination of nested integrals using Algorithm I, as described in Sec. 13 . 11 

Finally, since the upper bound of the outermost integral is one, the result can be 
expressed by MZVs with roots of unity (/3p = 1 with P t e N). 

Step 5 

The same as step 4 of Algorithm I (Sec. [2]). 



5 Applications of Algorithm II 

We apply Algorithm II for evaluating a 3-loop integral 

r d D P d D k d D K 1 

2 J (2ti) d (2ti) d (2tt) d (k-v){p-v)k 2 p 2 K 2 {k + q + K) 2 {p + q + k) 2 ' 1 J 



** We note that c\(t') and c^it') in eq. (|73| can be expressed by products of powers of t', (1 ± t'), 
(1 ± it'). For this reason, the prefactors can also be expressed by these factors. Thus, the absorption of 
the prefactor does not alter the form eq. (|57]l . 
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This is another master integral necessary for a computation of 03; the diagram is shown 
in Fig. [3] The notations are the same as in Sec. 13.31 

Using Mellin-Barnes integral representation, one can express M2 by K with five 
different arguments with (L,M,N) = (1,1,2). After expanding in e, six T functions 
remain uncanceled in the summand, for each K. We need to compute the expansion co- 
efficients of K around three different arguments (ao,i, 60,1; c o,i> ^o,i> e o,i> e o,2, /o,i, /o,2) = 
(i, 1; |, 1; |, 1, 1, 2), (i, 1; 1, |; 1, §, f , §), (1, |; 1, |; |, 2, 2, 3) up to order e 5 , e 3 , e 3 , respec- 
tively. We can use Algorithm II to evaluate these coefficients and obtain 



M 2 = A(e)(-g 2 )- 3e x 



28vr 4 f 226C(5) 116vr 2 C(3) 4 / , 224 \ 

J +tt 4 4log2 1 



135e I 3 9 V 135/ 



3 3 V V2/ 9 3 

, 1792 20 log 2 2 , \ 428tt 6 

~\~ 7T + -^ + 32l0g2 + 



135 3 J 2835 



, f / ~\ \ j- / \ 480 w , 9l 384 , 384 3072 

e 2 |768Li 4 (i)C(3) + — C(3) 2 log2 + 1536s 6 - — s 6 log2 + — s 7a + —s 7b 

4960 C (7) 14464 C (5) + 64^ + 32C(3) ^ g _ ^ g 



7T 



21 3 3 

2 ' 



+ 



(-5121^) + 128Li 5 (|) - - S _ 32M2 

M 



16 log 5 2 _ 64 4 \ 7r4 /31457C(3) 14336 40 log 3 2 
15 T ° g / + V 945 135 - 

160 log 2 2 



(60) 



where s 7a = z(-5, 1, l)+z(-6, l)+z(-5, 2)+z(-7) = £ m >„> fc >i ^ = -0.9529600757562 

and s 7b = z(5, -1, -l)+z(7)+z(5, 2)+z(-6, -1) = E m >„> fc >i = 1.029121262964324 . 

We have checked the above result numerically by evaluating Feynman parameter integral 
representation using sector decomposition. This result has not been published elsewhere. 

We have computed yet another master integral for S3 using Algorithm II. This is 
I\q given in [32] Jj where the analytical expressions for the coefficients of the e _1 and 
e° terms are presented. Ref. [32] uses the Mellin-Barnes method in combination with 
PSLQ algorithm [33] (a sophisticated estimate of analytical results from high-precision 
numerical results) for the computation. We have reproduced their result, thus providing 
a cross check in a fully analytical way. 



* Ref. [35] presents analytical results for the expansion coefficients of (selected) 16 master integrals 
for (X3. 
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6 Summary and Discussion 



We have constructed algorithms for computing two types of multiple sums (Algorithms 
I and II) which appear in higher-order loop computations. For instance, these types of 
sums appear in analytic evaluations of Mellin-Barnes integral representations by closing 
contours and expanding in e. 

Algorithm I applies to a multiple sum without F functions [eq. ©]: 



bN(u,h,...,%N-i) 



bi(v) 62(^,11) 

sm= E E ■ 

il=a\{v) i2=a 2 (v,ii) i N =a N (u,i 1 ,...,i N _ 1 ) 



E 



\1N 

N 



l\ r L r (u,ix, 



in) 



(61) 



The algorithm is valid for arbitrary multiplicity N and works efficiently (This is partly 
due to our specific regularization scheme for treating surface terms at infinity; see dis- 
cussion at the end of Sec. 2.) The results are expressed by Z-sums [with argument nv 
(n G N)] and generalized MZVs. 

There are benefits for constructing an efficient algorithm specifically for this type 
of sums. First, this type of sums appear frequently in computation of expansion co- 
efficients of loop integrals in e. Secondly, it has other useful applications as shown in 
Sec. [3j (i) Conversion of non- nested integrals to nested integrals (used in Algorithm II), 
(ii) Deriving non-trivial relations among MZVs; for example, we find a further reduction 
of the proposed basis elements at (w, I) = (2, 8) in [T5] . 

Algorithm II is used to evaluate the expansion coefficients of a double sum with T 
functions [eq. (pH]) ]. 



K(a,b;c,d;e,f) = 



£,m=0 



n 



r (i+a p ) 



' M 

n 



r (m + c q ) 
r (m+d ) 



n x r(£+m+e r 
t T(£+m + f r 



(62) 



around half-integer values (including integer values) of its arguments (a ,b ; c ,d ; e ,f ). 
For p, q, r > 2, ao, p — &o,p, co, g — d 0tq , eo, r — fo,r are assumed to be integers. Hence, there 
remain at most six uncanceled T functions in the summand in each expansion coefficient. 
The algorithm applies to arbitrary expansion coefficients. The results are expressed by 
generalized MZVs with roots of unity. 

These algorithms are useful in evaluating some complicated master integrals which 
appear in computation of the 3-loop correction to the static QCD potential (a 3 ). We 
have presented new results for two master integrals (Sees. 13.31 and IB}. 

In practical implementation of the algorithm^, one can make various improvements 
to compute efficiently. For Algorithm I, however, at the moment we do not have ideas for 
essential modifications of the algorithm, apart from technical refinements. On the other 
hand, for Algorithm II, we may improve efficiencies non-trivially by categorizing the 
arguments of K and changing algorithms according to the categories. Although in most 



' A Mathematica package " Wo," for computing multiple sums using the algorithms developed in this 
paper is available at http://www.tuhep.phys.tohoku.ac. jp/^program/ with examples and instruc- 
tions. In particular, Type I sums can be evaluated fairly efficiently. 
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difficult cases K can be evaluated only by the algorithm described in Sec. HI in simpler 
cases, it is possible to construct algorithms which work substantially faster: for instance, 
we may express K in an integral representation involving hypergeometric functions 2 -Pi 
and utilize various identities of 2-^1- We explain details of our practical implementation 
elsewhere. 

In computation of 03 we encounter multiple sums which cannot be evaluated with 
Algorithms I and II. An example is a sum, where the arguments of V functions in the 
numerator and denominator include the summation indices in different combinations, 
such as: 

ET(m + a)T(n + b)T(k + m + c)T(k + m + n + d) ■ ■ ■ 
T(k + a')T(m + n + b')T(k + 2m + n + d) ■ ■ ■ ^ ' 

k,m,n,--- v ' 

To our knowledge, there are no known technologies which can be used to evaluate these 
sums in general cases. We are currently studying how to generalize the algorithms in 
order to evaluate these different types of sums. It may be useful to combine the present 
algorithms with other methods for analytical evaluation of loop integrals, such as the 
method of differential equation, glue-and-cut method, etc. 



Appendices 

A Definitions and Conventions 

Z-sums are defined by 

Z (w; an, «2, • • • , a n ; Pi, P2, ■ ■ ■ , Pn) = 

■w>k\>k2>—>k n 



>x kTkT---k^ 



W nkl nk 2 fcr.-1-l nk n 

V^V — ••■ V — (64) 



Z-^t h ai £-~< h a2 t-^ h 

fcl=l 1 k 2 = l 2 fc«=l 



where a iy w G N and pi G C. For w = 00, the sums are (generalized) multiple zeta 
valued (MZVs) or multiple polylogarithms of Goncharov: 

Z (00; ai, a 2 , • • • , a n ; fix, p 2 , ■ ■ ■ , Pn) = Li an) ..., aa (Pn, ■ ■ ■ , Pi) ■ (65) 

We also use a short-hand notation for MZVs, in contexts where it is possible to restore 
the original forms: 

Z(oo; ai,..., a n ; Pi,..., P n ) = z(aiPi, a n p n ). (66) 



* It is more common to restrict to the case ft = ±1 when referring to MZVs. In view of applications 
to higher-order loop computations, where this restriction becomes inconvenient, we relax conditions on 
j3i in this paper. 
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Harmonic polylogarithms (HPLs) are defined recursively: For weight one, they are 
defined b}0 

/o (*) = -, f±i{x) = -!—, (67) 
x x =F 1 

#o (x) = log (x) , H ±1 (x)= [ dt f ±1 (t) = log(l=Fx), (68) 

Jo 

and for higher weights, 

1 /"* 

#o,-,o(a0 = -r [log(x)] n , H aua2i ^ an (x) = dt f ai (t)H(t;a 2 ,...,a n ), (69) 
v ^v-' ra! Jo 

n 

where at least one of ax, a 2 , ■ ■ ■ , a n is non-zero. The above definition can be generalized 
to the case 

f p (x) = — ; peC (70) 
x — p 

where a; of -ff a i,a 2 ,...,a n (x) becomes a general complex number. We refer to it as a general- 
ized HPL. H ai ,a 2 ,...,a„(x) is another representation of an MZV or multiple polylogarithm 
[eq. dS5}]i 



B Procedure for Raising or Decreasing Exponents 



In this appendix we explain the procedure used in step 4 of Algorithm II (see Sec. H]). 
It is used to raise or decrease the exponents of the factors in the integrand, until all of 
them become non-negative when the regularization parameters are set to zero. 
We consider an integral of the form 



dt 



P(log<p 1 (t,t'),...,log ( j> 9 (t,t'))H(t) 



(71) 



Here, 4> s (t,t') denotes one of the nine factors t, (1 ± £), (1 ± it), (1 ± t/t'), (1 ± 
tt'); P(xx, ■ ■ ■ ,xg) represents a polynomial of x\, . . . ,xg; H(t) is a generalized HPL of 
eqs. (!67|) - (!70|) ; 6 = (5x, ■ ■ ■ , S n ) denote the regularization parameters; the exponent a s (6) 
is a linear function of 8, where a s (0) is an integer. 

A rough sketch of the procedure is as follows. We first multiply the integrand by one, 
which can be decomposed as 



1 = t + (1 - t) = - + ( 1 - - 1 = (1 - t) + — ^— ( 1 + - 1 

v ; f v tJ 1 + t i + t'\ t'J 



(72) 



^ The definition in this paper differs from the conventional one, /±i(ir) = j^, considering the 
generalization eq. (1701) 

* The summation representation is obtained by first writing the integrand of nested integral repre- 
sentation in Taylor series and then integrating. 
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etc. This is repeated until the exponents are raised sufficiently. Then we integrate by 
parts to trade between different exponents. 

Let us explain in detail. First we select a pair of the factors <p s and <f) a , for which 
a s (0) and a a (0) are both negative. We multiply the integrand of eq. flTTj) by 

l = c s (t')<p s (t,t') + c a (t l )<p a (t,t'). (73) 

This is repeated |a s (0) + 0:0.(0) | times. After expanding the integrand, in each term, the 
exponents of (fi s and (fi a are both non-negative or at least one of them is non-negative|§ 
We repeat this process successively for pairs of negative exponents in each term. In the 
end, in each term, at most one negative exponent remains. At the same time, the sum 
of all the exponents are non-negative for each term. Next we perform integration by 
parts in every term which has <\> v s with v < 0, in such a way to integrate This is 



repeated until the exponent of cj) s becomes non-negative. In the end, all the exponents 
can be made non-negative. 

In the case of a double integral with respect to £1, £2, we iterate the above procedure 
twice, first with respect to the inner integral and then the outer integral. Due to the 
absence of the factors (l±£i/£ 2 ), (l±$i6) ; iteration of integration by parts with respect 
to £2 can render all the exponents non- negative. (The same is true for 771,772-) 
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